from . import query
from . import ifilter


import itk

# import gc
# 
# gc.disable()

def subtractImages(x,y):
    try:
        iType = query.getImageType(x)
        sub = itk.SubtractImageFilter[iType,iType,iType].New()
        sub.SetInput1(x)
        sub.SetInput2(y)
        sub.Update()
        #sub.thisown = 0
        img = sub.GetOutput()
        return img
    except Exception, error:
        print "failed in subtractImages", error
def divideImages(x,y):
    """divide the two images x and y. Must have the same voxel type and dimension"""
    try:
        iType = query.getImageType(x)
        mf = itk.DivideImageFilter[iType,iType,iType].New()
        mf.SetInput1(x)
        mf.SetInput2(y)
        mf.Update()
        #mf.thisown = 0
        
        img = mf.GetOutput()
        img.Update()
        return img
    except Exception, error:
        print "failed in multiplyImages", error
def multiplyImages(x,y):
    """multiply the two images x and y. Images must have the same voxel type and dimension"""
    try:
        xtype = query.getImageType(x)
        ytype = query.getImageType(y)
        if( ytype != xtype ):
            y2 = ifilter.castImage(y,oitype=xtype)
        else:
            y2 = y
        
            mf = itk.MultiplyImageFilter[xtype,xtype,xtype].New()
            mf.SetInput1(x)
            mf.SetInput2(y2)
            mf.Update()
            #mf.thisown = 0
            img = mf.GetOutput()
            return img
    except Exception, error:
        print "failed in multiplyImages", error
def sumImages(images):
    """Sum a series of images contained in images"""
    try:
        img = images[0]
        iType = query.getImageType(img)
        af = itk.AddImageFilter[iType,iType,iType].New()
        #af.thisown = 0
        for i in images[1:]:
            af.SetInput1(img)
            af.SetInput2(i)
            af.Update()
            img = af.GetOutput()   
        img.Update() 
        return img
    except Exception, error:
        print "failed in sumImages", error
def meanImage(images):
    """Compute the mean image for the series of images in images"""
    try:
        itype = query.getImageType(images[0])
        sumx = sumImages(images,iType=iType)
        img = scaleOffsetImage(sumx,1.0/len(images),0,iType=iType)
        return img
    except Exception, error:
        print "failed in meanImage", error
def createMatchedUniform(x,off,iType=itk.Image[itk.F,3], **args):
    try:
        sc = itk.ShiftScaleImageFilter[iType,iType].New()
        #sc.thisown = 0
        sc.SetInput(x)
        sc.SetShift(0)
        sc.SetScale(0) 
        #sc.Update()
        #temp = sc.GetOutput()

        sc2 = itk.ShiftScaleImageFilter[iType,iType].New()
        sc2.SetInput(sc.GetOutput())
        sc2.SetShift(off)
        sc2.SetScale(1.0)
        sc2.Update()
        #sc2.thisown = 0
        img = sc2.GetOutput() 
        return img    
    except Exception, error:
        print "failed in createMatchedUniform", error
def scaleOffsetImage(x,scale,off,iType=itk.Image[itk.F,3]):
    try:
        sc = itk.ShiftScaleImageFilter[iType,iType].New()
        sc.SetInput(x)
        sc.SetShift(off)
        sc.SetScale(scale)      
        sc.Update()
        #sc.thisown = 0
        img = sc.GetOutput()  
        return img
    except Exception, error:
        print "failed in scaleOffsetImage", error
def regressImages(y,xvals,iType = itk.Image[itk.F,3]):
    try:
        shape = query.getImgSize(y[0],returnITK=False)
        x = []

        for xv in xvals:
            print xv
            i = createMatchedUniform(y[0],xv)
            x.append(i)
        #print "sum images"
        sumx = sumImages(x,iType=iType)
        sumy = sumImages(y,iType=iType)
        
        xx = []
        xy = []
        for i in range(len(x)):
            #print "xx"
            xx.append(multiplyImages(x[i],x[i]))
            #print "yy"
            xy.append(multiplyImages(x[i],y[i]))
        sumxx = sumImages(xx,iType=iType)
        sumxy = sumImages(xy,iType=iType)
        sxsy = multiplyImages(sumx,sumy,iType=iType)
        sxsx = multiplyImages(sumx,sumx,iType=iType)

        meanX = scaleOffsetImage(sumx,1.0/len(x),0,iType=iType)
        mNum = subtractImages(sumxy,scaleOffsetImage(sxsy,1.0/len(x),0,iType=iType))
        mDen = subtractImages(sumxx,scaleOffsetImage(sxsx,1.0/len(x),0,iType=iType))
        
        bx = divideImages(mNum,mDen,iType=iType)
        ax = subtractImages(meanX,multiplyImages(bx,meanX,iType=iType))
        return bx,ax
    except Exception, error:
        print "failed in regressImages", error
        
            
        
        
        
    
